/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "reaction/reaction.h"
#include "molecule/molecule_arom_match.h"

Reaction::Reaction()
{
   clear();
}

void Reaction::clear() {
   _reactantCount = 0;
   _productCount = 0;
   _catalystCount = 0;
   _allMolecules.clear();
   _allFlags.clear();
   _indexes.clear();
   name.clear();
}

int Reaction::addReactantMolecule() {
   ++_reactantCount;
   return _addMolecule(REACTANT, 0, 0, 0);
}

int Reaction::addProductMolecule() {
   ++_productCount;
   return _addMolecule(PRODUCT, 0, 0, 0);
}

int Reaction::addCatalystMolecule() {
   ++_catalystCount;
   return _addMolecule(CATALYST, 0, 0, 0);
}

int Reaction::addReactantMoleculeCopy(const Molecule& mol, Array<int>* mapping, Array<int> *inv_mapping) {
   ++_reactantCount;
   return _addMolecule(REACTANT, &mol, mapping, inv_mapping);
}

int Reaction::addProductMoleculeCopy(const Molecule& mol, Array<int>* mapping, Array<int> *inv_mapping) {
   ++_productCount;
   return _addMolecule(PRODUCT, &mol, mapping, inv_mapping);
}

int Reaction::addCatalystMoleculeCopy(const Molecule& mol, Array<int>* mapping, Array<int> *inv_mapping) {
   ++_catalystCount;
   return _addMolecule(CATALYST, &mol, mapping, inv_mapping);
}

int Reaction::findMoleculeSideIdx(const Molecule* qmol) const
{
   int mol_idx = findMoleculeIdx(qmol);
   int side_idx = 0;
   int type = _indexes[mol_idx];

   for (int i = begin(); i < mol_idx; i = next(i))
      if(_indexes[i] == type)
         side_idx++;
   return side_idx;
}

int Reaction::findMoleculeIdx(const Molecule* qmol) const {
   for(int i = begin(); i < end(); i = next(i)) {
      if(qmol == &_allMolecules.at(i))
         return i;
   }
   throw Error("cannot find aam number");
}

int Reaction::findAtomByAAM (int mol_idx, int aam) const
{
   const Molecule &mol = _allMolecules.at(mol_idx);

   for (int i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
      if (getAAM(mol_idx, i) == aam)
         return i;
   return -1;
}

int Reaction::findAamNumber(const Molecule* qmol, int atom_number) const {
   for(int i = begin(); i < end(); i = next(i)) {
      if(qmol == &_allMolecules.at(i))
         return getAAM(i, atom_number);
   }
   throw Error("cannot find aam number");
}

int Reaction::findInversionNumber(const Molecule* qmol, int atom_number) const {
   for(int i = begin(); i < end(); i = next(i)) {
      if(qmol == &_allMolecules.at(i))
         return getInversion(i, atom_number);
   }
   throw Error("cannot find inversion");
}

int Reaction::findReactingCenter(const Molecule* qmol, int bond_number) const {
   for(int i = begin(); i < end(); i = next(i)) {
      if(qmol == &_allMolecules.at(i))
         return getReactingCenter(i, bond_number);
   }
   throw Error("cannot find reacting center");
}

int Reaction::findExactChange(const Molecule* qmol, int atom_number) const {
   for(int i = begin(); i < end(); i = next(i)) {
      if(qmol == &_allMolecules.at(i))
         return getExactChange(i, atom_number);
   }
   throw Error("cannot find exact change");
}

void Reaction::aromatizeBonds() {
   for (int i = begin(); i < end(); i = next(i)) {
      MoleculeAromatizer::aromatizeBonds(_allMolecules[i]);
   }
}

void Reaction::dearomatizeBonds() {
   for (int i = begin(); i < end(); i = next(i)) {
      MoleculeDearomatizer::dearomatizeMolecule(_allMolecules[i]);
   }
}

void Reaction::aromatizeQueryBonds() {
   for (int i = begin(); i < end(); i = next(i)) {
      QueryMoleculeAromatizer::aromatizeBonds(_allMolecules[i]);
   }
}

void Reaction::calcImplicitHydrogens(bool allow_bad) {
   for (int i = begin(); i < end(); i = next(i)) {
      _allMolecules[i].calcImplicitHydrogens(allow_bad);
   }
}

void Reaction::calcConnectivity() {
   for (int i = begin(); i < end(); i = next(i)) {
      _allMolecules[i].calcConnectivity(true);
   }
}

void Reaction::findBondsInRings() {
   for (int i = begin(); i < end(); i = next(i)) {
      _allMolecules[i].findBondsInRings();
   }
}

bool Reaction::isAllConnected() const {
   for (int i = begin(); i < end(); i = next(i)) {
      if (!Graph::isConnected(_allMolecules[i]))
         return false;
   }
   return true;
}

void Reaction::markStereocenterBonds() {
   for (int i = begin(); i < end(); i = next(i)) {
      _allMolecules[i].getStereocenters().markBonds();
   }
}

bool Reaction::haveCoord (const Reaction &reaction) {
   for (int i = reaction.begin(); i < reaction.end(); i = reaction.next(i)) {
      if (!reaction.getMolecule(i).have_xyz)
         return false;
   }
   return true;
}

void Reaction::makeTransposedForSubstructure(const Reaction &other) {

   clear();

   QS_DEF(Array<int>, transposition);
   int index = 0;

   for (int i = other.begin(); i < other.end(); i = other.next(i)) {

      //make transposition for molecules
      other._transposeMoleculeForSubstructure(i, transposition);
      
      switch(other._indexes[i]) {
         case REACTANT:
            index = addReactantMoleculeCopy(other.getMolecule(i), 0, 0);
            break;
         case PRODUCT:
            index = addProductMoleculeCopy(other.getMolecule(i), 0, 0);
            break;
         case CATALYST:
            index = addCatalystMoleculeCopy(other.getMolecule(i), 0, 0);
            break;
      }

      getMolecule(index).makeSubmolecule(other.getMolecule(i), transposition, 0);

      for (int j = 0; j < transposition.size(); j++) {
         getAAMArray(index).at(j) = other.getAAM(i, transposition[j]);
         getInversionArray(index).at(j) = other.getInversion(i, transposition[j]);
         getExactChangeArray(index).at(j) = other.getExactChange(i, transposition[j]);
      }

      for (int j = getMolecule(index).edgeBegin(); j < getMolecule(index).edgeEnd(); j = getMolecule(index).edgeNext(j)) {
         const Edge &edge = getMolecule(index).getEdge(j);
         int edge_idx = other.getMolecule(i).findEdgeIndex(transposition[edge.beg], transposition[edge.end]);
         getReactingCenterArray(index).at(j) = other.getReactingCenter(i, edge_idx);
      }
   }

}

void Reaction::cloneReaction(const Reaction& other, ObjArray< Array<int> >* mappings, ObjArray< Array<int> >* inv_mappings) {
   clear();
   int index = 0;

   QS_DEF(ObjArray< Array<int> >, tmp_mappings);
   ObjArray< Array<int> >* mol_mappings = mappings;

   if(mol_mappings == 0) {
      mol_mappings = &tmp_mappings;
   }
   mol_mappings->clear();
   if(inv_mappings != 0)
      inv_mappings->clear();

   for(int i = 0; i < other.end(); ++i) {
      mol_mappings->push();
      if(inv_mappings != 0)
         inv_mappings->push();
   }

   for(int i = other.begin(); i < other.end(); i = other.next(i)) {
      const Molecule& rmol = other.getMolecule(i);
      Array<int> * inv_mapping = 0;
      if(inv_mappings != 0)
         inv_mapping = &inv_mappings->at(i);
      switch(other._indexes[i]) {
         case REACTANT:
            index = addReactantMoleculeCopy(rmol, &mol_mappings->at(i), inv_mapping);
            break;
         case PRODUCT:
            index = addProductMoleculeCopy(rmol, &mol_mappings->at(i), inv_mapping);
            break;
         case CATALYST:
            index = addCatalystMoleculeCopy(rmol, &mol_mappings->at(i), inv_mapping);
            break;
      }

      for(int j = rmol.vertexBegin(); j < rmol.vertexEnd(); j = rmol.vertexNext(j)) {
         getAAMArray(index).at(j) = other.getAAM(i, mol_mappings->at(i)[j]);
         getInversionArray(index).at(j) = other.getInversion(i, mol_mappings->at(i)[j]);
      }
      //for query
      getExactChangeArray(index).resize(other.getExactChangeArray(i).size());
      if(getExactChangeArray(index).size() > 0) {
         for(int j = rmol.vertexBegin(); j < rmol.vertexEnd(); j = rmol.vertexNext(j)) {
            getExactChangeArray(index).at(j) = other.getExactChange(i, mol_mappings->at(i)[j]);
         }
      }

      for (int j = getMolecule(index).edgeBegin(); j < getMolecule(index).edgeEnd(); j = getMolecule(index).edgeNext(j)) {
         const Edge &edge = getMolecule(index).getEdge(j);
         int edge_idx = other.getMolecule(i).findEdgeIndex(mol_mappings->at(i)[edge.beg], mol_mappings->at(i)[edge.end]);
         getReactingCenterArray(index).at(j) = other.getReactingCenter(i, edge_idx);
      }

   }
}

int Reaction::_addMolecule(int type, const Molecule* mol, Array<int>* mapping, Array<int> *inv_mapping) {
   Molecule& new_mol = _allMolecules.push();
   if (mol != 0)
      new_mol.cloneMolecule(*mol, mapping, inv_mapping);
   
   _allFlags.push().init(new_mol);
   _indexes.push(type);

   return (_indexes.size()-1);
}

int Reaction::_nextElement(int type, int index) const {
   for (++index; index < _indexes.size(); ++index) {
      if (_indexes[index] & type)
         break;
   }
   return index;
}
                 
int Reaction::sideCount (int side) const { 
   switch (side) {
   case REACTANT:
      return reactantsCount();
   case PRODUCT: 
      return productsCount();
   case CATALYST:
      return catalystCount();
   default:
      throw Error("Unknown side identifier");
   }
}

int Reaction::_compare(const int &i1, const int &i2, const void *c) {
   const _SortingContext &context = *(const _SortingContext *) c;

   bool is_pseudo1 = context.m.isPseudoAtom(i1);
   bool is_pseudo2 = context.m.isPseudoAtom(i2);
   if (is_pseudo1 != is_pseudo2)
   {
      if (is_pseudo1)
         return -1;
      return 1;
   }
   
   // Compare by AAM, reacting centers and other reacting flags
   int res = context.rdata[i2] - context.rdata[i1];

   if (res != 0 || is_pseudo1)
      return res;

   // Compare by atom frequency
   int labels_by_freq[] = {ELEM_C, ELEM_H, ELEM_O, ELEM_N, ELEM_P, ELEM_F,
      ELEM_S, ELEM_Si, ELEM_Cl, ELEM_Br, ELEM_I, ELEM_At};

   int label1 = context.m.getAtom(i1).label;
   int label2 = context.m.getAtom(i2).label;
   int idx1, idx2;

   for (idx1 = 0; idx1 < (int) NELEM(labels_by_freq); idx1++)
      if (label1 == labels_by_freq[idx1])
         break;
   for (idx2 = 0; idx2 < (int) NELEM(labels_by_freq); idx2++)
      if (label2 == labels_by_freq[idx2])
         break;

   res = idx2 - idx1;

   if (res != 0)
      return res;

   // compare by degree

   return context.m.getVertex(i2).degree() - context.m.getVertex(i1).degree();
}

void Reaction::_transposeMoleculeForSubstructure(int index, Array<int> &transposition) const {
   QS_DEF(Array<int>, has_reacting_info);
   const Molecule& mol = _allMolecules[index];
   const _ReactionFlags& react_flags = _allFlags[index];

   has_reacting_info.clear_resize(mol.vertexEnd());
   has_reacting_info.zerofill();
   transposition.clear();

   for (int i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i)) {
      if (react_flags.atomAtomMapping.at(i) > 0)
         has_reacting_info[i] += 4;
      if (react_flags.inversionNumbers.at(i) > 0 || react_flags.queryExactChanges.at(i) > 0)
         has_reacting_info[i] += 1;

      transposition.push(i);
   }

   for (int i = mol.edgeBegin(); i < mol.edgeEnd(); i = mol.edgeNext(i)) {
      if (react_flags.reactingCenters.at(i) > 0) {
         const Edge &edge = mol.getEdge(i);
         has_reacting_info[edge.beg] += 2;
         has_reacting_info[edge.end] += 2;
      }
   }

   _SortingContext context(mol, has_reacting_info);

   transposition.qsort(_compare, &context);
}

void Reaction::clearAAM() {
   for (int i = begin(); i < end(); i = next(i)) {
      _allFlags[i].atomAtomMapping.zerofill();
   }
}

void ReactionHighlighting::init(const Reaction &reaction) {
   _graphHighlightings.clear();
   for(int i = 0; i < reaction.end(); i++) {
      _graphHighlightings.push();
   }
   for(int i = reaction.begin();i < reaction.end(); i = reaction.next(i))
      _graphHighlightings[i].init(reaction.getMolecule(i));
}

void ReactionHighlighting::clear() {
   _graphHighlightings.clear();
}

void ReactionHighlighting::nondestructiveInit(const Reaction &reaction) {
   while (_graphHighlightings.size() < reaction.end())
      _graphHighlightings.push();
}

void Reaction::_ReactionFlags::init(int vertex_count, int edge_count) {
   atomAtomMapping.resize(vertex_count);
   reactingCenters.resize(edge_count);
   inversionNumbers.resize(vertex_count);
   queryExactChanges.resize(vertex_count);

   atomAtomMapping.zerofill();
   reactingCenters.zerofill();
   inversionNumbers.zerofill();
   queryExactChanges.zerofill();
}

void ReactionAtomNeighbourhoodCounters::calculate(const Reaction &reac, bool is_query) {
   int i;

   _counters.resize(reac.count());

   for (i = reac.begin(); i < reac.end(); i = reac.next(i))
      _counters[i].calculate(reac.getMolecule(i), is_query);
}

void Reaction::saveBondTypes(const Reaction& reaction, ObjArray< Array<int> > &bond_types) {

   while (bond_types.size() < reaction.end())
      bond_types.push();

   int i;

   for (i = reaction.begin(); i != reaction.end(); i = reaction.next(i))  {
      Molecule::saveBondTypes(reaction.getMolecule(i), bond_types[i]);
   }
}

void Reaction::loadBondTypes(Reaction& reaction, ObjArray< Array<int> > &bond_types) {

   int i;

   for (i = reaction.begin(); i != reaction.end(); i = reaction.next(i))  {
      Molecule::loadBondTypes(reaction.getMolecule(i), bond_types[i]);
   }
}
